Structural stability of finite dispersion-relation 

preserving schemes 



Claire David^*, Pierre Sagaut^ 

Universite Pierre et Marie Curie-Paris 6 
^ Institut Jean Le Rond d'Alembert, UMR CNRS 7190, 
Boite courrier 71*^162, 4 place Jussieu, 75252 Paris, cedex 05, France 

April 18, 2008 



Abstract 

The goal of this work is to determine classes of travelling solitary wave so- 
lutions for a differential approximation of a finite difference scheme by means 
of a hyperbolic ansatz. It is shown that spurious solitary waves can occur in 
finite-difference solutions of nonlinear wave equation. The occurance of such 
a spurious solitary wave, which exhibits a very long life time, results in a non- 
vanishing numerical error for arbitrary time in unbounded numerical domain. 
Such a behavior is referred here to has a structural instability of the scheme, 
since the space of solutions spanned by the numerical scheme encompasses 
types of solutions (solitary waves in the present case) that are not solution 
of the original continuous equations. This paper extends our previous work 
about classical schemes to dispersion-relation preserving schemes. 



1 Introduction: The DRP scheme 

The Burgers equation: 

Ut + CUUx- ^lUxx = ^-, (1) 

c, being real constants, plays a crucial role in the history of wave equations. It 
was named after its use by Burgers [2] for studying turbulence in 1939. 
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i, n denoting natural integers, a linear finite difference scheme for this equation 
can be written under the form: 

Y^ai^uT = (2) 

where: 

ur = u{lh,mT) (3) 

/ G {i — 1, i, i + 1}, m G {n — 1, n, n + 1}, j = 0, n^, n = 0, rif. The aim 
are real coefficients, which depend on the mesh size h, and the time step r. 
The Courant-Friedrichs-Lewy number [cfl) is defined as a = cr/h . 
A numerical scheme is specified by selecting appropriate values of the coefficients 
aim- Then, depending on them, one can obtain optimum schemes, for which the 
error will be minimal. 



m being a strictly positive integer, the first derivative |^ is approximated at the 
j^th ^q(Iq Qf ^i^Q spatial mesh by: 

(4) 

k=—m 

Following the method exposed by C. Tam and J. Webb in [5], the coefficients 
gammak are determined requiring the Fourier Transform of the finite difference 
scheme (jll) to be a close approximation of the partial derivative ( |^ )«• 
(jlj) is a special case of: 



^a^)' - Yl lku{x + kh) (5) 

k=—m 

where x is a continuous variable, and can be recovered setting x = I h. 

Denote by uj the phase. Applying the Fourier transform, referred to by ^, to both 

sides of ([5]), yields: 

m 

jujuc^ ^ -^ke^^'^^u (6) 

k=—m 

j denoting the complex square root of —1. 

Comparing the two sides of ([6]) enables us to identify the wavenumber A of the finite 

m 

difference scheme dl]) and the quantity 4 TfcC-'^'^^, i. e.: The wavenumber of 

k=—m 

the finite difference scheme (jl]) is thus: 
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,jkujh 



(7) 



k=—m 



To ensure that the Fourier transform of the finite difference scheme is a good 
approximation of the partial derivative over the range of waves with wave- 

length longer than 4 h, the a priori unknowns coefficients 7^ must be choosen so as 
to minimize the integrated error: 

£ = J^^\Xh-Xh\'^d{Xh) 

^ m 

k=—m 

m 

= IL K + j E {co«( ^ C) + J sin( k 0} I' dC 



k=-r 



C- Y I'k sin( k C) 

k=—m 
m 

C- Y Ik sin( k C) 



k——m 

The conditions that £ is a minimum are: 
i. e.: 



E 7fe cos( fc C) 

_k——m 

m 

E 7fe cos( /c C) 



.k=—m 



dC 
dC 



-Csin(iC) + XI cos((A;-i)C) > = 

k=—m ) 

Changing i into — i, and k into — /c in the summation yields: 

C sin(iC) + 1-k cos ( (-/c + i) C) WC = 



fe=— m 



1. e.: 



C sin(iC) + J] 7-fc cos ( (A; - i) C) > ciC = 



k=—m 



Thus: 



(8) 



(9) 



(10) 



(12) 
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/i — 'II' 

I 5Z {7-fe + 7fc} cos((A;-OC) rfC = (13) 

•^0 k=-m 

which yields: 

kj^i,k=—m ^ 

which can be considered as a hnear system of 2 m + 1 equations, the unknowns of 
which are the 7_j + 7j, i = —m, . . . , m. The determinant of this system is not 
equal to zero, while it is the case of its second member: the Cramer formulae give 
then, for i = —m, . . . , m: 

7-^ + 7^ = (15) 

or: 



For i = 0, one of course obtains: 



All this ensures: 



7-i = -7i (16) 



70 = (17) 



k=—m 

The values of the •jk coefficients are obtained by substituting relations (fT6il into 
(HUD: 

m 

5^ 7fc = (19) 

k=—m 

m being a strictly positive integer, a 2m + 1-points DRP scheme ([5]) is thus given 
by: 



,n+l I I 



m 

r 



ur'^ + u:-+^ ^ 7fc<+fe = o (20) 



k— — r 



where the 7^, k G {— m,m} are the coefficients of the considered scheme, and 
satisfy the relations (fT6i) . 



4 



Considering again the ui"^ terms as functions of the mesh size h and time step r, 
expanding them at a given order by means of their Taylor series expansion, and 
neglecting the o(t^) and o{h'^) terms, for given values of the integers p, q, leads to 
the following differential approximation (see [H]): 

+ J E ^^^^'c(-,||,|^,/^,r)=0 (21) 

where J-'l^^j^ denotes the function of m, h, r obtained by means of the above 

Taylor expansion, r, s being integers. 

For sake of simplicity, a non-dimensional form of Eq. (12T1) will be used: 

- ir* + + E 7. (s. 0. ^) = (22) 

k=—m ^ 

Depending on this differential approximation fl22p . solutions, as solitary waves, may 
arise. 

The paper is organized as follows. DRP schemes are analyzed in section [21 The 
general method is exposed in Section [31 Classical DRP schemes are studied in 
section [H where it is shown that out of the two studied schemes, only one leads to 
solitary waves. A related class of travelling wave solutions of equation f[?T|) is thus 
presented, by using a hyperbolic ansatz. 

2 Analysis of DRP schemes 

Consider function of the time step r, and expand it at the second order 

by means of its Taylor series: 

Ui"~^^ = u {i h, (n + l) t) = u {i h,nT) +t Ut {i h,nT) + — Uu {i h,nT) + o(r^) (23) 
It ensures: 

— — = Ut{ih,nT) + - utt {i h.nr) + o{t) (24) 

T I 

In the same way, for k G {— m, m}, consider Mj+fc'^ as a function of the mesh size 
/i, and expand it at the fourth order by means of its Taylor series expansion: 
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iti-)-fe" =u {{i + k) h, n t) 

= u {h,n t) + k hux {i h,nT) . . 

2 2 4 4 
H — Uxx (i h,nT) A u^xxx (i /i, n t) + o(h'*) 



Equation ( l2Ti) can thus be written as: 



■ ut {i h,n t) — — utt {i /i, n r) + o(r) 



^ ™ r ^2^2 ^4^4 ^ 
+ - 22 Ik Y^{ih,nT) + khux {ih,nT) -\ — Uxx{ih,nT)A Uxxxx {i h,nT) + o{h*) j = 



(26) 



k— — m 

i. e., at X = i h and t = nr: 



- Ut - ^ Utt + o{t) + ^ -Tklu + khUxA hUxx — Uxxxx + o{h'^)\ = (27) 

2 Ai , ^ — ' I 2 4! J 

T9|) ensures then: 

m 

- Mt - -Utt +o(t) + - ^ fc7fc {/lu^ +o(/i'')} = (28) 



fc— — m 



The related first differential approximation of the Burgers equation ([T]) is thus 
obtained neglecting the o(r) and o(/i^) terms, yielding: 

m 

- Mt - - Mtt + r ^ k^i-Ux =0 (29) 

fc— — m 

For sake of simplicity, this latter equation can be adimensionalized in the following 

way: 

set: 

u = UqU 

t = Tot (30) 
where: 

f/o = ^ (31) 

To 

In the following, Rch will denotes the mesh Reynolds number, defined as: 

Unh , , 

Rch = ^ 32 
/i 



For j G IN, the change of variables (!30|) leads to: 
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Ut 
Up 



(33) 



becomes: 



Uo T Uo y—^ ^0 

— 'H - 77 :;2 "tt + 2 T 2^ A; 7fc — «s = (34) 



Multiplying (El) by ^ yields: 



T T TQ yr-^ 

ui - - — + 2 -— y k-ykhu^ = (35) 
2, ff-' i~> '■' 



For h = ho, due to cr = Eq. (l35l) becomes: 



r/io 



'.TO tJ.Reh 

which simplifies in: 



T T Ho V — r 

+ = (36) 



fj 2(7 

■ - o "tt + '^k'ykUi=0 (37) 



3 Solitary waves 

Approximated solutions of the Burgers equation ([T]) by means of the difference 
scheme ( l20l) strongly depend on the values of the time and space steps. For specific 
values of r and h, equation flHTl) can, for instance, exhibit travelling wave solutions 
which can represent great disturbances of the searched solution. 
We presently aim at determining the conditions, depending on the values of the 
parameters r and h, which give birth to travelling wave solutions of ( 1371) . 
Following Feng [3] and our previous work in which travelling wave solutions 
of the CBKDV equation were exhibited as combinations of bell- profile waves and 
kink-profile waves, we aim at determining travelling wave solutions of fl37p (see 

Following [3], we assume that equation (1371) has travelling wave solutions of the 
form 

?i(x,t)=M(0, i = x-vi (38) 
where v is the wave velocity. Substituting (IHSjl into equation (1^ leads to: 

^{u,u^'\{-vyu^'^) = Q, (39) 
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Performing an integration of ( l39i) with respect to C, leads to an equation of the 
form: 

= C, (40) 

where C is an arbitrary integration constant, which will be the starting point for 
the determination of solitary waves solutions. 

It is important to note that, contrary to other works, the integration constant is 
not taken equal to zero, which would lead to a loss of solutions. 

4 Travelling Solitary Waves 

4.1 Hyperbolic Ansatz 

The discussion in the preceding section provides us useful information when we 
construct travelling solitary wave solutions for equation (IH^ . Based on these re- 
sults, in this section, a class of travelling wave solutions of the equivalent equation 
( l29l) is searched as a combination of bell-profile waves and kink-profile waves of the 
form 

n 

u{x, i)=J2 {U^ tanh^ [C,{x-vi)]+ Vi sech' [a{x - vi + xq)]) + Vq (41) 
1=1 

where the U^s, V-s, C-s, {i = I, ■ ■ ■ , n), Vq and v are constants to be determined. 
In the following, c is taken equal to 1. 

4.2 Theoretical analysis 

Substitution of (13^ into equation (HUj) leads to an equation of the form 

J2 A tanh* {Q ^) sech^' {d ^) sinh^ {d ^) = C (42) 

the Ai being real constants. 

The difficulty for solving equation (H2l) lies in finding the values of the constants Ui, 
Vi, Ci, Vq and v by solving the over-determined algebraic equations. Following [3], 
after balancing the higher-order derivative term and the leading nonlinear term, 
we deduce n = 1. 
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Then, following [4j we replace sech(Ci^) by 7crq~=crT; sinh(Ci^) by — — f — —, 

tanh(Ci ^) by gc|g^g-cio multiply both sides by (l + e^^*"^)^, so that equation 
fj42|) can be rewritten in the following form: 

4 

J^PkiUu Vu Ci, V, \/o)e'=^^« = 0, (43) 

fc=0 

where the (/c = 0, 4), are polynomials of Ui, Vi, Ci, Vq and v. 

Depending wether ( H2|) admits or no consistent solutions, spurious solitary waves 

solutions may, or not, appear. 



4.3 Numerical scheme analysis 

Equation (l39l) is presently given by: 



+ ^fe7fe«'(0 = (44) 

2 Mfie^ ^1 



Performing an integration of (jSj) with respect to yields: 



1. e.: 



■^m-^n'{0 + ^^J2k^ku{^) = C (45) 



^fc^,. _„U(0- ^ = C (46) 



I 2 



where C is an arbitrary integration constant. 
Substitution of fl52|) for = 1 into equation P6|) leads to: 



2(T 



1. e.: 



{{/i tanh [ Ci 5] + soch [ Ci 5] + Vq} - — 



l/isech2[Ci5]- Vi 



sinhfCig] 
cosh^ [Cxi\ 



■■ C 



(47) 



2<7 ™ If e^'iS-e-Ci? 2Vi 1 i''' f / 2 eC'i^-e-Ci? 1 



l^^Jeh j I eCi«+e-Ci« eCi« + e-Ci« ' "J 2 [ ^ ^ e"^! « + e""^! V (e Ci ? + e - Ci 5) 

2 ^^^^ 

Multiplying both sides by (l + e^*"^^) yields: 
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which is a fourth-order equation in e^^^. This equation being satisfied for any real 
value of ^, one therefore deduces that the coefficients of e'^'"^^, A; = 0, . . . , 4 must 
be equal to zero, i.e.: 



2cr 



^A:7fc-«[> {-Ui + Vo}-"^-^ {4Ui + 2Vi} = C 



2(7 
2(7 
2(7 



J^k-lk - v \ 2Vi=0 



J^i^yk-vi Vo = o 

k-1 J 
m '\ 

^kik-v \ Vi + Ci (T Vi = 



(50) 



2(7 

IxRch 



^fe7fe-«MC^i + ^o} = 



V — 



2(7 

IJ-R^h 



/c 7fc , Vi 7^ leads to the trivial nul solution. Therefore, Vi is neces- 



k=i 



sarily equal to zero, which implies: 



2 (7 



IxReu 
Ui = - 



k=l 

c 



(51) 



2 Ci 1)2 (T 

Vo en , Ci en 



All Z)i?P schemes admit thus kink-profile travelling solitary waves solutions, given 
by: 



u{x,i) = — - 



C 



■tanh[Ci(:r-t;i)] -|- 



(52) 



5 Conclusions 

The analysis of the nonlinear equivalent differential equation for finite-differenced 
DRP schemes for the Burgers equation has been carried out. We show that all DRP 
schemes admit spurious travelling solitary waves solutions, which make them, as 
regards this point, structurally instable. 
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